1 Introduction

This is a clean analysis of the R code in here:

/Volumes/Groupdir/SUN-DAN-Brickman/Jose/Elena_clean/RNAseq_analyses/nascent_RNAseq_20200526/scripts/DESEQ2_analysis.R

We will make a DESEQ2 analysis focusing only on the 8h timepoint

2 Load libraries

library(DESeq2)
library(tidyverse)
library(RColorBrewer) #colour scheme for hierarchical clustering heatmap
library(ggrepel) #labels for ggplot
library(ggdendro)
library(reshape2)
library(apeglm) #for MAplot function in DESEq2
library(genefilter) #top var genes
library(pheatmap)
library(readxl)
library(org.Mm.eg.db)

Create results folder for this notebook

results_folder <- "../results/DESEQ2_analysis/"

if(!file.exists(results_folder)){
  dir.create(results_folder)
}

3 Useful functions

3.1 PCA biplot from rlog transformation

my_PCAbiplot <- function(object, ...) {
  .local <- function (object,  X = "PC1", Y = "PC2", intgroup = "condition", ntop = 500, nplot = 10, 
                      returnData = FALSE) {
    
    rv <- rowVars(assay(object))
    select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, 
                                                       length(rv)))]
    genes <- mcols(object)$GeneSymbol[select]
    
    pca <- prcomp(t(assay(object)[select, ]))
    row.names(pca$rotation) <- genes
    
    percentVar <- pca$sdev^2/sum(pca$sdev^2)
    
    if (!all(intgroup %in% names(colData(object)))) {
      stop("the argument 'intgroup' should specify columns of colData(dds)")
    }
    
    intgroup.df <- as.data.frame(colData(object)[, intgroup, 
                                                 drop = FALSE])
    group <- if (length(intgroup) > 1) {
      factor(apply(intgroup.df, 1, paste, collapse = ":"))
    }
    else {
      colData(object)[[intgroup]]
    }
    d <- data.frame(pca$x, group = group, 
                    intgroup.df, name = colnames(object))
    
    if (returnData) {
      attr(d, "percentVar") <- percentVar
      return(d)
    }
    
    plot <- ggplot(data = d, aes_string(x = X, y = Y, color = "group")) + 
      geom_point(size = 3) + xlab(paste0(X,": ", round(percentVar[which(colnames(pca$x)==X)] * 100), "% variance")) + 
      ylab(paste0("Y",": ", round(percentVar[which(colnames(pca$x)==Y)] * 100), "% variance")) + theme_bw()
    
    plot <- plot + geom_hline(yintercept = 0, size=.2) + geom_vline(xintercept = 0, size=.2)
    
    datapc <- data.frame(varnames=rownames(pca$rotation), pca$rotation)
    mult <- min(
      (max(d[,Y]) - min(d[,Y])/(max(datapc[,Y])-min(datapc[,Y]))),
      (max(d[,X]) - min(d[,X])/(max(datapc[,X])-min(datapc[,X])))
    )
    datapc <- transform(datapc,
                        v1 = .7 * mult * (get(X)),
                        v2 = .7 * mult * (get(Y))
    )
    datapc$distance <- sqrt(datapc$v1^2 + datapc$v2^2)
    #datapc <- datapc[select,]
    
    plot <- plot + coord_equal() + geom_text_repel(data=datapc[1:nplot,], aes(x=v1, y=v2, label=varnames), size = 3, vjust=1, color="black")
    plot <- plot + geom_segment(data=datapc[1:nplot,], aes(x=0, y=0, xend=v1, yend=v2), arrow=arrow(length=unit(0.2,"cm")), alpha=0.75, color="blue")
    plot
  }
  .local(object, ...)
}

3.2 MA and volcano plot with LFC shrinkage

DE_plot_func <- function(dds, contrast ,LFC, sign, out) {
  
  res_LFC <- lfcShrink(dds, coef = (contrast), type = "apeglm")
  
  res_LFC <- as.data.frame(res_LFC)
  res_LFC$gene_name <- mcols(dds)$GeneSymbol
  res_LFC[["padj"]][is.na(res_LFC[["padj"]])]<-0.99
  res_LFC[["log2FoldChange"]][is.na(res_LFC[["log2FoldChange"]])]<-0
  res_LFC$sig <- "No"
  res_LFC$sig[abs(res_LFC$log2FoldChange) > LFC & res_LFC$padj < sign] <- "Yes"
  res_LFC$sig <- factor(res_LFC$sig, levels <- c("Yes","No"))
  
  MA_plot<-ggplot(res_LFC,aes(x=baseMean, y=log2FoldChange,col=sig))+
    geom_point(alpha=0.5,size=0.8)+
    geom_hline(yintercept=0,alpha=0.75,col="red")+
    labs(col="")+
    xlab("Base Mean")+
    ylab("Log2 Fold Change") +
    scale_color_manual(labels=c("DE","Not DE"),values=c("red","black"))+ xlim(0,1000)
  theme_bw()
  
  volcano_plot <- ggplot(res_LFC, aes(x=log2FoldChange, y = -log10(padj))) +
    geom_jitter(aes(color = sig)) + theme_bw() + 
    geom_hline(yintercept=-log10(sign), linetype="dashed", color = "black") + 
    geom_vline (xintercept = c(-LFC, LFC), linetype = "dashed", color = "black") + 
    guides(color = FALSE) + ylab("-log10(Adjusted p-value)") + xlab("Log2 FC")
  
  if (out == FALSE){
    return(list(MA_plot, volcano_plot))
    
  }else{
    ggsave(plot = MA_plot, filename = paste0(out,"_MAplot.pdf"), 
         height = 7, width = 7)
    
    ggsave(plot = volcano_plot, filename = paste0(out,"_volcano.pdf"), 
         height = 7, width = 7)
  }
}

4 Load data

data_path <- ("../data/counts_gene-id_gene.txt")
data <- read.delim(data_path,as.is=T)

Annotations of the data. “-11” is to pick the 8h samples (which are the 12 last samples of the file)

annot <- data[,1:6]

counts <- as.matrix(data[,(ncol(data)-11):ncol(data)])

Getting gene names from the mouse ensembl data base

genes <- annot$Geneid
db <- org.Mm.eg.db
ensembl <- suppressWarnings(mapIds(db, keys=genes, keytype="ENSEMBL", column="SYMBOL",multiVals="first"))
genes <- data.frame(SYMBOL = ensembl, ENSEMBL = genes, stringsAsFactors = F)
genes$SYMBOL[is.na(genes$SYMBOL)] = genes$ENSEMBL[is.na(genes$SYMBOL)]

annot$GeneSymbol = genes$SYMBOL

Getting sample annotations for DEseq analysis. We extract information (metadata) from sample column names

names = colnames(counts)
names = sapply(names,function(x) gsub("X8h","X_8h",x))
p1 = "(.*)_(.*)_(.*)_(.*)_(.*)_uniq"
colData <- data.frame(str_match(names,p1),stringsAsFactors = F)
colnames(colData) = c("Sample","Dox","Time","Tam","Rep","Rep2")
colData$SampleName = colnames(counts)
colData$TimePt = as.numeric(str_match(colData$Time,"^(.*)h")[,2])

Samples to fix: the Tam field should be empty and the values there should move to Dox Jose’s NOTE: This does not seem to apply to the 8h time point, i believe it is for the ones before the 2nd treatment

idx = colData$Dox=="X"
colData[idx,"Dox"] = colData[idx,"Tam"]
colData[idx,"Tam"] = NA

colData$SampleNameShort = with(colData,paste(Dox,Tam,Time,Rep2, sep="_"))
colData$ShortLabel = with(colData,paste(Dox,Tam,Time, sep="_"))
colData$ShortLabelTime = with(colData,paste(Dox,Tam, sep="_"))

5 Differential expression using H2O_ETOH as reference

colData$condition = factor(colData$ShortLabel, levels = c("H2O_ETOH_8h","H2O_4OHT_8h","DOX_ETOH_8h","DOX_4OHT_8h"))
colData$Dox = factor(colData$Dox, levels = c("H2O","DOX"))
colData$Tam = factor(colData$Tam, levels = c("ETOH", "4OHT"))

5.1 DESeq analysis

#full condition
dds <- DESeqDataSetFromMatrix(counts, colData, design = ~ condition)

5.2 Pre-filtering the dataset

First we add the extracted annotation that we did at the beginning and add it to the information that we just retrieved

mcols(dds) = cbind(annot,mcols(dds))
nrow(dds)
## [1] 25783
dds <- dds[rowSums(counts(dds)) > 10, ]
nrow(dds)
## [1] 17798
dds <- DESeq(dds)

5.3 Exploratory analysis

5.3.1 Rlog transformation

rld <- rlog(dds)

rv <- rowVars(assay(rld))
select <- order(rv, decreasing = TRUE)[seq_len(min(500, length(rv)))]
genes <- mcols(rld)$GeneSymbol[select]

5.3.2 PCA plot

pca <- prcomp(t(assay(rld)[select, ]))
row.names(pca$rotation) <- genes

percentVar <- pca$sdev^2/sum(pca$sdev^2)
d <- data.frame(pca$x, pErk = as.character(colData$Tam), Oct4 = as.character(colData$Dox),stringsAsFactors = F)
d$pErk[d$pErk == "4OHT"] <- "+"
d$pErk[d$pErk == "ETOH"] <- "-"
d$Oct4[d$Oct4 == "H2O"] <- "+"
d$Oct4[d$Oct4 == "DOX"] <- "-"
PCA_plot <- ggplot(data = d, aes_string(x = "PC1", y = "PC2", color = "Oct4", shape = "pErk")) + 
  geom_point(size = 3) + xlab(paste0("PC1",": ", round(percentVar[which(colnames(pca$x)=="PC1")] * 100), "% variance")) + 
  ylab(paste0("PC2",": ", round(percentVar[which(colnames(pca$x)=="PC2")] * 100), "% variance")) + theme_bw()

PCA_plot <- PCA_plot + geom_hline(yintercept = 0, size=.2) + geom_vline(xintercept = 0, size=.2) + coord_fixed(ratio=1.5) 

PCA_plot

PCA plot with loadings. Can make plots of other PC. nplot plots the n genes with most variance. We are not using it here

#my_PCAbiplot(rld, X = "PC1", Y = "PC2", nplot = 20) + ggtitle("PCA on rlg transformation and loadings")

5.3.3 Sample distances from rlog transformation

sampleDists <- dist(t(assay(rld)))
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- rld$Rep2
colnames(sampleDistMatrix) <- rld$Rep2
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
annot_col <- colData[,c("Rep2", "Dox","Tam")]
rownames(annot_col) <- annot_col$Rep2
annot_col$pErk[annot_col$Tam == "4OHT"] <- "+"
annot_col$pErk[annot_col$Tam == "ETOH"] <- "-"
annot_col$Oct4[annot_col$Dox == "H2O"] <- "+"
annot_col$Oct4[annot_col$Dox == "DOX"] <- "-"
annot_col <- annot_col[,-c(1,2,3)]
ann_colors <- list(Oct4 = c("+"="grey45", "-"= "grey75"),pErk = c("+"="grey45", "-"= "grey75"))
sample_distances <- pheatmap(sampleDistMatrix, show_rownames = F, show_colnames = F,
         clustering_distance_rows=sampleDists,
         clustering_distance_cols=sampleDists,
         col=colors, annotation_col = annot_col, annotation_colors = ann_colors)

6 DE analysis 4OHT vs Control

##MA and volcano plot with LFC shrinkage

plots_4OHT_vs_Control <- DE_plot_func(dds, "condition_H2O_4OHT_8h_vs_H2O_ETOH_8h", 1, 0.05, out = FALSE)
plots_4OHT_vs_Control[1]
## [[1]]
plots_4OHT_vs_Control[2]
## [[1]]

6.1 Up and down regulated pErk regulated genes

res <- results(dds, name = "condition_H2O_4OHT_8h_vs_H2O_ETOH_8h")
res <- as.data.frame(res)
res$gene_name <- make.names(mcols(dds)$GeneSymbol, unique = T)
res[["padj"]][is.na(res[["padj"]])]<-0.99
res[["log2FoldChange"]][is.na(res[["log2FoldChange"]])]<-0

up<-filter(res,log2FoldChange>1, padj<0.05)
down<-filter(res,log2FoldChange<(-1),padj<0.05)

tam_regulated <- filter(res, abs(log2FoldChange)> 1, padj < 0.05)
tam_regulated$DE <- "Up"
tam_regulated$DE[tam_regulated$log2FoldChange < (-1)] <- "Down"

Write results of DE genes

write.table(tam_regulated, file = paste0(results_folder,"4OHT_vs_Control_DE_genes.tsv"), sep = "\t",
            quote = F, col.names = T, row.names = F)

6.2 Heatmap of genes regulated by pErk

#mat <- assay(normTransform(dds))
mat <- assay(rlog(dds))
rownames(mat) <- res$gene_name
mat  <- mat[tam_regulated$gene_name,]
mat <- mat[,c("H2O_8h_ETOH_a_S19_uniq", "H2O_8h_ETOH_b_S20_uniq", "H2O_8h_ETOH_d_S21_uniq",
              "DOX_8h_ETOH_a_S22_uniq", "DOX_8h_ETOH_b_S23_uniq", "DOX_8h_ETOH_d_S24_uniq",
              "H2O_8h_4OHT_a_S25_uniq", "H2O_8h_4OHT_b_S26_uniq", "H2O_8h_4OHT_d_S27_uniq",
              "DOX_8h_4OHT_a_S28_uniq", "DOX_8h_4OHT_b_S29_uniq", "DOX_8h_4OHT_d_S30_uniq")]

Row annotation (gene annotation)

row_anno <- data.frame(row.names = rownames(mat), pErk_DE = tam_regulated$DE)
anno <- as.data.frame(colData(dds)[, c("Dox","Tam")])

colnames(anno) <- c("Oct4", "pErk")
anno$Oct4 <- as.character(anno$Oct4)
anno$Oct4[which(anno$Oct4 == "DOX")] <- "-"
anno$Oct4[which(anno$Oct4 == "H2O")] <- "+"
anno$Oct4 <- as.factor(anno$Oct4)

anno$pErk <- as.character(anno$pErk)
anno$pErk[which(anno$pErk == "4OHT")] <- "+"
anno$pErk[which(anno$pErk == "ETOH")] <- "-"
anno$pErk <- as.factor(anno$pErk)

ann_colors <- list(Oct4 = c("+"="grey45", "-"= "grey75"),pErk = c("+"="grey45", "-"= "grey75"),
                   pErk_DE = c("Up" = "firebrick2", "Down" = "skyblue2"))
hm <- pheatmap(mat, annotation_col = anno, annotation_colors = ann_colors, 
              annotation_row = row_anno,
              show_rownames = F, show_colnames = F, 
              cluster_cols = F, cluster_rows = T, 
              main = "Dysregulation of Fgf-Erk dependent genes\nAbsolute LFC >1, adj. p-value <0.05", scale = "row",
              treeheight_row = 0, treeheight_col = 0)

7 DE analysis DOX vs Control

7.1 MA and volcano plot with LFC shrinkage

plots_DOX_vs_Control <- DE_plot_func(dds, "condition_DOX_ETOH_8h_vs_H2O_ETOH_8h", 1, 0.05, out = FALSE)
plots_DOX_vs_Control[1]
## [[1]]
plots_DOX_vs_Control[2]
## [[1]]

7.2 Up and down regulated genes

res <- results(dds, name = "condition_DOX_ETOH_8h_vs_H2O_ETOH_8h")
res <- as.data.frame(res)
res$gene_name <- mcols(dds)$GeneSymbol
res[["padj"]][is.na(res[["padj"]])]<-0.99
res[["log2FoldChange"]][is.na(res[["log2FoldChange"]])]<-0

up<-filter(res,log2FoldChange>1, padj<0.05)
down<-filter(res,log2FoldChange<(-1),padj<0.05)

dox_regulated <- filter(res, abs(log2FoldChange)> 1, padj < 0.05)
dox_regulated$DE <- "Up"
dox_regulated$DE[dox_regulated$log2FoldChange < (-1)] <- "Down"

Write results of DE genes

write.table(dox_regulated, file = paste0(results_folder,"DOX_vs_Control_DE_genes.tsv"), sep = "\t",
            quote = F, col.names = T, row.names = F)

8 DE analysis DOX-4OHT vs Control

8.1 MA and volcano plot with LFC shrinkage

plots_DOX_4OHT_vs_Control <- DE_plot_func(dds, "condition_DOX_4OHT_8h_vs_H2O_ETOH_8h", 1, 0.05, out = FALSE)
plots_DOX_4OHT_vs_Control[1]
## [[1]]
plots_DOX_4OHT_vs_Control[2]
## [[1]]

8.2 Up and down regulated genes

res <- results(dds, name = "condition_DOX_4OHT_8h_vs_H2O_ETOH_8h")
res <- as.data.frame(res)
res$gene_name <- mcols(dds)$GeneSymbol
res[["padj"]][is.na(res[["padj"]])]<-0.99
res[["log2FoldChange"]][is.na(res[["log2FoldChange"]])]<-0

up<-filter(res,log2FoldChange>1, padj<0.05)
down<-filter(res,log2FoldChange<(-1),padj<0.05)

tam_dox_regulated <- filter(res, abs(log2FoldChange)> 1, padj < 0.05)
tam_dox_regulated$DE <- "Up"
tam_dox_regulated$DE[tam_dox_regulated$log2FoldChange < (-1)] <- "Down"

Write results

write.table(dox_regulated, file = paste0(results_folder, "DOX_4OHT_vs_Control_DE_genes.tsv"), sep = "\t",
            quote = F, col.names = T, row.names = F)

9 Data frame for venn diagram

tam_dox_up <- tam_dox_regulated$gene_name[tam_dox_regulated$DE == "Up"]
dox_up <- dox_regulated$gene_name[dox_regulated$DE == "Up"]
tam_up <- tam_regulated$gene_name[tam_regulated$DE == "Up"]
maximum <- max(length(tam_dox_up), length(dox_up), length(tam_up))
up <- data.frame(DOX_TAM = c(tam_dox_up,rep("",maximum - length(tam_dox_up))),
                 DOX = c(dox_up,rep("",maximum - length(dox_up))),
                 TAM = c(tam_up,rep("",maximum - length(tam_up)))
                 )
write.table(up, file = paste0(results_folder, "chow_ruskey_UP_data.tsv"), sep = "\t",
            quote = F, col.names = T, row.names = F)
tam_dox_down <- tam_dox_regulated$gene_name[tam_dox_regulated$DE == "Down"]
dox_down <- dox_regulated$gene_name[dox_regulated$DE == "Down"]
tam_down <- tam_regulated$gene_name[tam_regulated$DE == "Down"]
maximum <- max(length(tam_dox_down), length(dox_down), length(tam_down))
down <- data.frame(DOX_TAM = c(tam_dox_down,rep("",maximum - length(tam_dox_down))),
                 DOX = c(dox_down,rep("",maximum - length(dox_down))),
                 TAM = c(tam_down,rep("",maximum - length(tam_down)))
)
write.table(down, file = paste0(results_folder, "chow_ruskey_DOWN_data.tsv"), sep = "\t",
            quote = F, col.names = T, row.names = F)

10 DE analysis DOX_4OHT_vs_H2O_4OHT

colData relevel

colData$condition = factor(colData$ShortLabel, levels = c("H2O_4OHT_8h","H2O_ETOH_8h","DOX_ETOH_8h","DOX_4OHT_8h"))
colData$Dox = factor(colData$Dox, levels = c("H2O","DOX"))
colData$Tam = factor(colData$Tam, levels = c("ETOH", "4OHT"))
#full condition
dds <- DESeqDataSetFromMatrix(counts, colData, design = ~ condition)

#Adding the extracted annotation that we did at the beginning and add it to the information that we just retrieved
mcols(dds) = cbind(annot,mcols(dds))

## Pre-filtering the dataset
nrow(dds)
## [1] 25783
dds <- dds[rowSums(counts(dds)) > 10, ]
nrow(dds)
## [1] 17798
dds <- DESeq(dds)

10.1 MA and volcano plot after LFC shrinkage

plots_DOX_4OHT_vs_H2O_4OHT <- DE_plot_func(dds, "condition_DOX_4OHT_8h_vs_H2O_4OHT_8h", 1, 0.05, out = FALSE)
plots_DOX_4OHT_vs_H2O_4OHT[1]
## [[1]]
plots_DOX_4OHT_vs_H2O_4OHT[2]
## [[1]]

10.2 Up and down regulated genes

res <- results(dds, name = "condition_DOX_4OHT_8h_vs_H2O_4OHT_8h")
res <- as.data.frame(res)
res$gene_name <- mcols(dds)$GeneSymbol
res[["padj"]][is.na(res[["padj"]])]<-0.99
res[["log2FoldChange"]][is.na(res[["log2FoldChange"]])]<-0

up<-filter(res,log2FoldChange>1, padj<0.05)
down<-filter(res,log2FoldChange<(-1),padj<0.05)

up_and_down_dox_tam_vs_h2o_tam <- filter(res, abs(log2FoldChange)> 1, padj < 0.05)
up_and_down_dox_tam_vs_h2o_tam$DE <- "Up"
up_and_down_dox_tam_vs_h2o_tam$DE[up_and_down_dox_tam_vs_h2o_tam$log2FoldChange < (-1)] <- "Down"

#gene_order <- order(up_and_down_dox_tam_vs_h2o_tam$log2FoldChange, decreasing = T)
#up_and_down_dox_tam_vs_h2o_tam <- up_and_down_dox_tam_vs_h2o_tam$gene_name[gene_order]

Write results

write.table(up_and_down_dox_tam_vs_h2o_tam, file = paste0(results_folder, "DOX_4OHT_vs_H2O_4OHT_DE_genes.tsv"), sep = "\t",
            quote = F, col.names = T, row.names = F)

11 Heatmap with clustering

## Gene clustering
tam_dox_and_tam <- intersect(up_and_down_dox_tam_vs_h2o_tam$gene_name, tam_regulated$gene_name)
#tam_dox_and_tam_only <- setdiff(tam_dox_and_tam, dox_regulated)
#tam_dox_and_tam_and_dox <- intersect(tam_dox_and_tam, dox_regulated)

genes <-mcols(rld)$GeneSymbol %in% tam_dox_and_tam
mat  <- assay(rld)[genes,]
mat <- mat[,c("H2O_8h_ETOH_a_S19_uniq", "H2O_8h_ETOH_b_S20_uniq", "H2O_8h_ETOH_d_S21_uniq",
              "DOX_8h_ETOH_a_S22_uniq", "DOX_8h_ETOH_b_S23_uniq", "DOX_8h_ETOH_d_S24_uniq",
              "H2O_8h_4OHT_a_S25_uniq", "H2O_8h_4OHT_b_S26_uniq", "H2O_8h_4OHT_d_S27_uniq",
              "DOX_8h_4OHT_a_S28_uniq", "DOX_8h_4OHT_b_S29_uniq", "DOX_8h_4OHT_d_S30_uniq")]
anno <- as.data.frame(colData(rld)[, c("Dox","Tam")])

colnames(anno) <- c("Oct4", "pErk")
anno$Oct4 <- as.character(anno$Oct4)
anno$Oct4[which(anno$Oct4 == "DOX")] <- "-"
anno$Oct4[which(anno$Oct4 == "H2O")] <- "+"
anno$Oct4 <- as.factor(anno$Oct4)


anno$pErk <- as.character(anno$pErk)
anno$pErk[which(anno$pErk == "4OHT")] <- "+"
anno$pErk[which(anno$pErk == "ETOH")] <- "-"
anno$pErk <- as.factor(anno$pErk)


rownames(mat) <- mcols(rld)$GeneSymbol[genes]

ann_colors <- list(Oct4 = c("+"="grey45", "-"= "grey75"),pErk = c("+"="grey45", "-"= "grey75"))
hm<- pheatmap(mat, annotation_col = anno, annotation_colors = ann_colors, show_rownames = F, show_colnames = F, 
              cluster_cols = F, cutree_row = 5, 
              main = "Dysregulation of Fgf-Erk dependent genes\nAbsolute LFC >1, adj. p-value <0.05", 
              cellwidth = 40)

11.1 Different clusters

clusters <- hclust(dist(mat, method = "euclidean"),method = "complete")

cluster_cut <- cutree(clusters, k = 5)

df <- data.frame (name = names(cluster_cut), cluster = cluster_cut)
write.table(df, file = paste0(results_folder,"genes_per_cluster.tsv"), sep = "\t", row.names = F, quote = F)
ncluster <- c(1:5)

for (i in ncluster){
  genes <- mcols(rld)$GeneSymbol %in% names(cluster_cut[which(cluster_cut == i)])
  mat  <- assay(rld)[genes,]
  
  mat <- mat[,c("H2O_8h_ETOH_a_S19_uniq", "H2O_8h_ETOH_b_S20_uniq", "H2O_8h_ETOH_d_S21_uniq",
                "DOX_8h_ETOH_a_S22_uniq", "DOX_8h_ETOH_b_S23_uniq", "DOX_8h_ETOH_d_S24_uniq",
                "H2O_8h_4OHT_a_S25_uniq", "H2O_8h_4OHT_b_S26_uniq", "H2O_8h_4OHT_d_S27_uniq",
                "DOX_8h_4OHT_a_S28_uniq", "DOX_8h_4OHT_b_S29_uniq", "DOX_8h_4OHT_d_S30_uniq")]
  anno <- as.data.frame(colData(rld)[, c("Dox","Tam")])
  rownames(mat) <- mcols(rld)$GeneSymbol[genes]
  pheatmap(mat, annotation_col = anno, cluster_cols = F, show_colnames = F, 
           main = paste0("Dysregulation of Fgf-Erk dependent genes. Gene cluster ", i,
                         "\nAbsolute LFC >1, adj. p-value <0.05 "), 
           cellwidth = 30)
}

12 Session info

devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.1.3 (2022-03-10)
##  os       macOS Monterey 12.4
##  system   x86_64, darwin21.3.0
##  ui       unknown
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       Europe/Copenhagen
##  date     2022-07-14
##  pandoc   2.18 @ /Applications/RStudio.app/Contents/MacOS/quarto/bin/tools/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package              * version    date (UTC) lib source
##  annotate               1.72.0     2021-10-26 [1] Bioconductor
##  AnnotationDbi        * 1.56.2     2021-11-09 [1] Bioconductor
##  apeglm               * 1.16.0     2021-10-26 [1] Bioconductor
##  assertthat             0.2.1      2019-03-21 [1] CRAN (R 4.1.2)
##  backports              1.4.1      2021-12-13 [1] CRAN (R 4.1.2)
##  bbmle                  1.0.24     2021-08-05 [1] CRAN (R 4.1.2)
##  bdsmatrix              1.3-4      2020-01-13 [1] CRAN (R 4.1.2)
##  Biobase              * 2.54.0     2021-10-26 [1] Bioconductor
##  BiocGenerics         * 0.40.0     2021-10-26 [1] Bioconductor
##  BiocParallel           1.28.3     2021-12-09 [1] Bioconductor
##  Biostrings             2.62.0     2021-10-26 [1] Bioconductor
##  bit                    4.0.4      2020-08-04 [1] CRAN (R 4.1.2)
##  bit64                  4.0.5      2020-08-30 [1] CRAN (R 4.1.2)
##  bitops                 1.0-7      2021-04-24 [1] CRAN (R 4.1.2)
##  blob                   1.2.2      2021-07-23 [1] CRAN (R 4.1.2)
##  brio                   1.1.3      2021-11-30 [1] CRAN (R 4.1.2)
##  broom                  0.7.12     2022-01-28 [1] CRAN (R 4.1.2)
##  bslib                  0.3.1      2021-10-06 [1] CRAN (R 4.1.2)
##  cachem                 1.0.6      2021-08-19 [1] CRAN (R 4.1.2)
##  callr                  3.7.0      2021-04-20 [1] CRAN (R 4.1.2)
##  cellranger             1.1.0      2016-07-27 [1] CRAN (R 4.1.2)
##  cli                    3.2.0      2022-02-14 [1] CRAN (R 4.1.2)
##  coda                   0.19-4     2020-09-30 [1] CRAN (R 4.1.2)
##  colorspace             2.0-3      2022-02-21 [1] CRAN (R 4.1.2)
##  crayon                 1.5.0      2022-02-14 [1] CRAN (R 4.1.2)
##  DBI                    1.1.2      2021-12-20 [1] CRAN (R 4.1.2)
##  dbplyr                 2.1.1      2021-04-06 [1] CRAN (R 4.1.2)
##  DelayedArray           0.20.0     2021-10-26 [1] Bioconductor
##  desc                   1.4.0      2021-09-28 [1] CRAN (R 4.1.2)
##  DESeq2               * 1.34.0     2021-10-26 [1] Bioconductor
##  devtools               2.4.3      2021-11-30 [1] CRAN (R 4.1.2)
##  digest                 0.6.29     2021-12-01 [1] CRAN (R 4.1.2)
##  dplyr                * 1.0.8      2022-02-08 [1] CRAN (R 4.1.2)
##  ellipsis               0.3.2      2021-04-29 [1] CRAN (R 4.1.2)
##  emdbook                1.3.12     2020-02-19 [1] CRAN (R 4.1.2)
##  evaluate               0.15       2022-02-18 [1] CRAN (R 4.1.2)
##  fansi                  1.0.2      2022-01-14 [1] CRAN (R 4.1.2)
##  farver                 2.1.0      2021-02-28 [1] CRAN (R 4.1.2)
##  fastmap                1.1.0      2021-01-25 [1] CRAN (R 4.1.2)
##  forcats              * 0.5.1      2021-01-27 [1] CRAN (R 4.1.2)
##  fs                     1.5.2      2021-12-08 [1] CRAN (R 4.1.2)
##  genefilter           * 1.76.0     2021-10-26 [1] Bioconductor
##  geneplotter            1.72.0     2021-10-26 [1] Bioconductor
##  generics               0.1.2      2022-01-31 [1] CRAN (R 4.1.2)
##  GenomeInfoDb         * 1.30.1     2022-01-30 [1] Bioconductor
##  GenomeInfoDbData       1.2.7      2022-03-03 [1] Bioconductor
##  GenomicRanges        * 1.46.1     2021-11-18 [1] Bioconductor
##  ggdendro             * 0.1.23     2022-02-16 [1] CRAN (R 4.1.2)
##  ggplot2              * 3.3.5      2021-06-25 [1] CRAN (R 4.1.2)
##  ggrepel              * 0.9.1      2021-01-15 [1] CRAN (R 4.1.2)
##  glue                   1.6.2      2022-02-24 [1] CRAN (R 4.1.2)
##  gtable                 0.3.0      2019-03-25 [1] CRAN (R 4.1.2)
##  haven                  2.4.3      2021-08-04 [1] CRAN (R 4.1.2)
##  highr                  0.9        2021-04-16 [1] CRAN (R 4.1.2)
##  hms                    1.1.1      2021-09-26 [1] CRAN (R 4.1.2)
##  htmltools              0.5.2      2021-08-25 [1] CRAN (R 4.1.2)
##  httr                   1.4.2      2020-07-20 [1] CRAN (R 4.1.2)
##  IRanges              * 2.28.0     2021-10-26 [1] Bioconductor
##  jquerylib              0.1.4      2021-04-26 [1] CRAN (R 4.1.2)
##  jsonlite               1.8.0      2022-02-22 [1] CRAN (R 4.1.2)
##  KEGGREST               1.34.0     2021-10-26 [1] Bioconductor
##  knitr                  1.37       2021-12-16 [1] CRAN (R 4.1.2)
##  labeling               0.4.2      2020-10-20 [1] CRAN (R 4.1.2)
##  lattice                0.20-45    2021-09-22 [2] CRAN (R 4.1.3)
##  lifecycle              1.0.1      2021-09-24 [1] CRAN (R 4.1.2)
##  locfit                 1.5-9.5    2022-03-03 [1] CRAN (R 4.1.2)
##  lubridate              1.8.0      2021-10-07 [1] CRAN (R 4.1.2)
##  magrittr               2.0.2      2022-01-26 [1] CRAN (R 4.1.2)
##  MASS                   7.3-55     2022-01-16 [2] CRAN (R 4.1.3)
##  Matrix                 1.4-0      2021-12-08 [2] CRAN (R 4.1.3)
##  MatrixGenerics       * 1.6.0      2021-10-26 [1] Bioconductor
##  matrixStats          * 0.61.0     2021-09-17 [1] CRAN (R 4.1.2)
##  memoise                2.0.1      2021-11-26 [1] CRAN (R 4.1.2)
##  modelr                 0.1.8      2020-05-19 [1] CRAN (R 4.1.2)
##  munsell                0.5.0      2018-06-12 [1] CRAN (R 4.1.2)
##  mvtnorm                1.1-3      2021-10-08 [1] CRAN (R 4.1.2)
##  numDeriv               2016.8-1.1 2019-06-06 [1] CRAN (R 4.1.2)
##  org.Mm.eg.db         * 3.14.0     2022-03-04 [1] Bioconductor
##  pheatmap             * 1.0.12     2019-01-04 [1] CRAN (R 4.1.2)
##  pillar                 1.7.0      2022-02-01 [1] CRAN (R 4.1.2)
##  pkgbuild               1.3.1      2021-12-20 [1] CRAN (R 4.1.2)
##  pkgconfig              2.0.3      2019-09-22 [1] CRAN (R 4.1.2)
##  pkgload                1.2.4      2021-11-30 [1] CRAN (R 4.1.2)
##  plyr                   1.8.6      2020-03-03 [1] CRAN (R 4.1.2)
##  png                    0.1-7      2013-12-03 [1] CRAN (R 4.1.2)
##  prettyunits            1.1.1      2020-01-24 [1] CRAN (R 4.1.2)
##  processx               3.5.2      2021-04-30 [1] CRAN (R 4.1.2)
##  ps                     1.6.0      2021-02-28 [1] CRAN (R 4.1.2)
##  purrr                * 0.3.4      2020-04-17 [1] CRAN (R 4.1.2)
##  R6                     2.5.1      2021-08-19 [1] CRAN (R 4.1.2)
##  RColorBrewer         * 1.1-2      2014-12-07 [1] CRAN (R 4.1.2)
##  Rcpp                   1.0.8      2022-01-13 [1] CRAN (R 4.1.2)
##  RCurl                  1.98-1.6   2022-02-08 [1] CRAN (R 4.1.2)
##  readr                * 2.1.2      2022-01-30 [1] CRAN (R 4.1.2)
##  readxl               * 1.3.1      2019-03-13 [1] CRAN (R 4.1.2)
##  remotes                2.4.2      2021-11-30 [1] CRAN (R 4.1.2)
##  reprex                 2.0.1      2021-08-05 [1] CRAN (R 4.1.2)
##  reshape2             * 1.4.4      2020-04-09 [1] CRAN (R 4.1.2)
##  rlang                  1.0.1      2022-02-03 [1] CRAN (R 4.1.2)
##  rmarkdown              2.12       2022-03-02 [1] CRAN (R 4.1.2)
##  rprojroot              2.0.2      2020-11-15 [1] CRAN (R 4.1.2)
##  RSQLite                2.2.10     2022-02-17 [1] CRAN (R 4.1.2)
##  rstudioapi             0.13       2020-11-12 [1] CRAN (R 4.1.2)
##  rvest                  1.0.2      2021-10-16 [1] CRAN (R 4.1.2)
##  S4Vectors            * 0.32.3     2021-11-21 [1] Bioconductor
##  sass                   0.4.0      2021-05-12 [1] CRAN (R 4.1.2)
##  scales                 1.1.1      2020-05-11 [1] CRAN (R 4.1.2)
##  sessioninfo            1.2.2      2021-12-06 [1] CRAN (R 4.1.2)
##  stringi                1.7.6      2021-11-29 [1] CRAN (R 4.1.2)
##  stringr              * 1.4.0      2019-02-10 [1] CRAN (R 4.1.2)
##  SummarizedExperiment * 1.24.0     2021-10-26 [1] Bioconductor
##  survival               3.2-13     2021-08-24 [2] CRAN (R 4.1.3)
##  testthat               3.1.2      2022-01-20 [1] CRAN (R 4.1.2)
##  tibble               * 3.1.6      2021-11-07 [1] CRAN (R 4.1.2)
##  tidyr                * 1.2.0      2022-02-01 [1] CRAN (R 4.1.2)
##  tidyselect             1.1.2      2022-02-21 [1] CRAN (R 4.1.2)
##  tidyverse            * 1.3.1      2021-04-15 [1] CRAN (R 4.1.2)
##  tzdb                   0.2.0      2021-10-27 [1] CRAN (R 4.1.2)
##  usethis                2.1.5      2021-12-09 [1] CRAN (R 4.1.2)
##  utf8                   1.2.2      2021-07-24 [1] CRAN (R 4.1.2)
##  vctrs                  0.3.8      2021-04-29 [1] CRAN (R 4.1.2)
##  withr                  2.5.0      2022-03-03 [1] CRAN (R 4.1.2)
##  xfun                   0.30       2022-03-02 [1] CRAN (R 4.1.2)
##  XML                    3.99-0.9   2022-02-24 [1] CRAN (R 4.1.2)
##  xml2                   1.3.3      2021-11-30 [1] CRAN (R 4.1.2)
##  xtable                 1.8-4      2019-04-21 [1] CRAN (R 4.1.2)
##  XVector                0.34.0     2021-10-26 [1] Bioconductor
##  yaml                   2.3.5      2022-02-21 [1] CRAN (R 4.1.2)
##  zlibbioc               1.40.0     2021-10-26 [1] Bioconductor
## 
##  [1] /usr/local/lib/R/4.1/site-library
##  [2] /usr/local/Cellar/r/4.1.3/lib/R/library
## 
## ──────────────────────────────────────────────────────────────────────────────